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One of the most fascinating experimental 
achievements of the last decade was the real- 
ization of Bose-Einstein Condensation (BEC) of 
ultra-cold atoms in optical lattices (OL's) [H [2j 
|3j |4]. The extraordinary level of control over 
these structures allows us to investigate complex 
solid state phenomena (H [Sj [U [3 HJ [9] and the 
emerging field of "atomtronics" promises a new 
generation of nanoscale devices. It is therefore 
of fundamental and technological importance to 
understand their dynamical properties. Here we 
study the outgoing atomic flux of BECs loaded in 
a one dimensional OL with leaking edges, using 
a mean fleld description provided by the Discrete 
Non-Linear Schrodinger Equation (DNLSE). We 
demonstrate that the atom population inside the 
OL decays in avalanches of size J. For intermedi- 
ate values of the interatomic interaction strength 
their distribution P{J) follows a power law i.e. 
7^(J) ~ 1/J" characterizing systems at phase tran- 
sition. This scale free behaviour of V{J) reflects 
the complexity and the hierarchical structure of 
the underlying classical mixed phase space. Our 
results are relevant in a variety of contexts (when- 
ever DNLSE is adequate), most prominently the 
light emmitance from coupled non-linear optics 
waveguides [12]. 

An optical lattice with a controlled leakage of the 
atomic BEC can be realized experimentally by the action 
of two separate continous microwave or Raman lasers to 
locally spin-flip BEC atoms (at the edges of the OL) to 
a nontrapped state jlOl [TT] . The spin-flipped atoms do 
not experience the magnetic trapping potential and are 
released through gravity in two atomic beams at the ends 
of the OL. The mathematical model (see Method section) 
that describes the dynamics of the BEC in a leaking OL 
of size M is 
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where 7 describes the atomic losses and = IV'nP 
is the atomic population at site n. Below we study 
the decay (due to leakage) of the total atom population 
■^ij) = X^n ^n{T) inside the OL as a function of the ini- 
tial effective interaction strength i.e. A — xA^(t = 0)/M. 

An exciting result appearing in the frame of nonlin- 
ear lattices is the existence of stationary, spatially local- 
ized solutions, termed Discrete Breathers (DB), which 
emerge due to the nonlinearity and discreteness of the 



system. DBs were observed in various experimental se- 
tups HH [H [171 [H [19] [20] while their existence and 
stability were studied thoroughly during the last decade 
[211 122I l23l l24l [25] . Their importance was already recog- 
nized in [26] where it was shown that they act as virtual 
bottlenecks which slow down the relaxation processes in 
generic nonlinear lattices |25l l26l |27] . Further works 
[27] established the fact that absorbing boundaries take 
generic initial conditions towards self-trapped DB's. Re- 
cently the same scenario was proposed for a BEC in a 
leaking OL where it was observed [TO] that N{t) decays 
in sudden bursts J (see inset of Fig. 1). Here, for the 
flrst time we present a full theoretical study of the decay 
process of N{t) and analyze the distribution V{J). We 
flnd three types of dynamical behaviour [ij: For very 
weak interactions A, the population decay is a smooth 
process which does not involve any avalanches. As inter- 
actions are increased, avalanches are created. For strong 
interaction strength, their distribution P( J) is exponen- 
tial. In contrast, for intermediate values of the interac- 
tion strength, P(J) follows a power law indicating the 
existence of a phase transition. Below, we will focus our 
analysis in this critical regime. 

Since we are interested in the effects of DBs on the re- 
laxation process, we introduce a localization parameter 
VTZ which provides a rough estimate of the relative num- 
ber of sites that are occupied by the remaining atoms in 
a leaking OL. It is deflned as 
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which in case of 7 = is the standard participation 
ratio. Accordingly, the more evenly the atoms spread 
over the lattice, the closer VTZ is to a constant of or- 
der 1. The VIZ approaches two limiting values p^4J: (a) 
VTZ = 1/2 corresponding to a random superposition of 
uncoupled sine waves for U = A — Q (linear regime) and 
(b) VTZ = 5/9 corresponding to a conflguration of atoms 
trapped in M uncoupled (T ^ ) wells for K ^ 00, 
which might be viewed as the formation of 0{M) DB's 
(multi-breather regime). In the closed system, the tran- 
sition between these limiting cases is smooth (see lower 
left inset of Fig. 1). 

In the open system (7 > 0) 7^7?, is a time dependent 
quantity. We numerically study the value 'PTis='PTZ{T*) 
for large times r* when the system has evolved from its 
initial thermalized state into a quasi-steady state [31] . 
For nonzero 7, instead of a smooth transition between 
the two extremes we observe a sharp drop of the VTZs at 
a critical interaction strength of A;, « 0.15 resembling a 
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phase transition (our numerics indicate that this transi- 
tion becomes a step function in the Hmit M <x>). VTZs 
drops down to its lowest possible value (c>c 1/M) corre- 
sponding to a single occupied site, i.e. the final state 
consists of a one self-trapped DB. If A increases further 
(Figs. 2b-c) the final state consists of a successively larger 
number of breathers, while for some critical A it evolves 
into a more complicated state with a power-law distribu- 
tion of norms Nn (see Fig. 1 left upper inset) |32] . 

In the following it is important to realize that if a 
breather solution exists for some value of A, it exists for 
all A's > (and for large enough M). This conclusion 
can be easily drawn by noting that a self-trapped DB is 
not directly coupled to the leaking edges, thus we can 
assume 7 = and then appropriately scale Eq. [5] There- 
fore breather solutions in particular do exist for A < Af, 
as well. So, what is the nature of the sharp transition 
at Afc? To answer this question, we have to examine the 
thermalized random initial states more closely. For small 
A w these states have exponentially distributed norms 
Pn^{x) = Mexp{—Mx) [11]. In addition, our numerical 
calculations (see the phase space analysis below) indi- 
cate that the minimal total nonlinearity associated with 
the central site of a self-trapped DB solution satisfy the 
relation MA\tfj\^ « 1 . Thus, a necessary condition in 
order to excite such a DB is that at least one site of the 
thermalized initial state has a norm x > xq = 1/(MA). 
Assuming independent random norms for the individ- 
ual sites we get that the probability for such an event 
is Wm^x) = 1 — (1 — exp{—Mx))^. By increasing M, 
we find that Wm{x) posseses a steep transition from 1 to 
at approximately X112 with VpM(a^i/2) = 1/2- We can 
thus give a rough estimation of the lower bound for Ab 
by demanding 2:1/2 > which leads to 
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in fair agreement with the value Af, w 0.15 for the system 
sizes we studied |33] . 

The most interesting situation emerges for large 
enough system sizes and (for the system sizes we could 
numerically study) interatomic interactions in the range 
of A « 0.5 — 1 j3l]. In this case we have found that 
the atoms leak out of the lattice in avalanches follow- 
ing a scale-free distribution, i.e. 'P{J) ~ while the 
norms Nn at individual lattice sites are distributed as 
Pn^{x) ~ x~'^with aw/? (see Fig. 1 left upper inset). 
This power law behavior for a whole range of parameter 
values A could be interpreted as a signature of a self- 
organized critical state |28j. 

Let us study in more detail the dynamics that lead to 
the creation of an avalanche in this parameter regime. 
One such event is depicted in Fig. 2d. A moving DB, 
coming from the bulk of the lattice, collides with the 
outer most self-trapped DB. As a result the self-trapped 
breather is shifted inwards by one site while at the same 



time a particle density proportional to the density of the 
moving breather tunnels through the self-trapped DB. 
Eventually this particle density will reach the leaking 
edge of the OL and decay in a form of an avalache. 

We have also found numerically [14] that during the 
migration process of the self-trapped DB, the number 
of particles and the energy of the three lattice sites in- 
volved is conserved, thus allowing us to turn the problem 
to the analysis of a reduced M — 3 system with inter- 
action strengths A in the critical range. As shown in 
Fig. 3a (inset), the non- linear trimer exhibits a hierar- 
chical mixed phase space structure with islands of regular 
motion (tori) embedded in a sea of chaotic trajectories. 
Trajectories inside the islands correspond to self-trapped 
DBs, provided that their frequency is outside the linear 
spectrum. In contrast, chaotic trajectories have contin- 
uous Fourier spectra, parts of which overlap with the 
linear spectrum of the infinite lattice ^1]. As long as the 
self-trapped DB is stable, it acts as a barrier which pre- 
vents atoms from reaching the leaking boundary. Thus, a 
necessary condition for an avalanche event is the destabi- 
lization of the DB. This is caused by a lattice excitation 
with particle density A''f'^'^*greater than a critical value 
which can push the regular orbit out of the island. At 
the same time a portion N^°'^ oc A^f '^'^'transmits through 
(see Fig. 3a). As can be seen from Fig. 2, the migra- 
tion process usually moves the DB inwards by one site 
which correspond to another island in the phase space of 
the extended lattice. Details of the migration process are 
still under investigation |141 130] . 

We conjecture that the size of avalanches J is propor- 
tional to the size S of islands in the mixed phase space 
found in the Poincare section of the non-linear trimer. A 
heuristic model [H] that mimics the hierarchical ('island- 
over island') structure of a typical mixed phase space 
leads us to a power law distribution P(S) ~ S~" where 
1 < a < 3. The lower bound comes from the requirement 
of having an infinite number of islands (self-similar prop- 
erty) while the upper-bound results from the requirement 
that the total volume of the phase space be finite. To ver- 
ify the above prediction for P{S) is a computationally 
demanding task. Therefore we used a numerically more 
convenient model: the kicked rotor which is a paradig- 
matic model of mixed phase space dynamics [29]. Leav- 
ing aside the technical details [H], we present in Fig. 3b 
the outcome of our numerical calculations. The results 
confirm the power law scaling of the island sizes over 
more than three orders of magnitude thus confirming the 
validity of our conjecture. 

In conclusion, we have shown that there is a critical 
regime of interatomic interactions, where particles are 
ejected out of a leaky OL in scale free avalanches. The 
observed power law distribution, is dictated by the hier- 
archical structure of the mixed phase space. Our results 
are quite universal and should be observable in many 
other experimental realizations of the DNLSE including 
molecular crystals, globular proteins and non-linear op- 
tics. 
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METHODS 

The simplest model that captures the dynamics of a 
dilute gas of bosonic atoms in a deep OL with chemical 
potential small compared to the vibrational level spac- 
ing, is the Bose-Hubbard Hamiltonian. In the case of 
weak interatomic interactions (superfiuid limit) and/or a 
large number of atoms per well (so that the total num- 
ber of atoms N ~ 0(10'' — 10^) is much bigger than the 
number of wells M), a further simplification is available 
since the BECs dynamics admits a semiclassical mean 
field description [13]. The resulting Hamiltonian is 

M M-1 

n = J2[u\iJ^\' + ^lM^]~-Y,irni'n+l + c.c.) (4) 

n— 1 n— 1 

where n is the well index, |?/;„(t)p = Nn{t) is the mean 
number of bosons at well n, U = iTrfi^asVes/rn de- 
scribes the interaction between two atoms on a single 
site {VcB is the effective mode volume of each site, m is 
the atomic mass, and is the s-wave atomic scattering 
length), /i„ is the well chemical potential, and T is the 
tunneling amplitude. The "wavefunction amplitudes" 
tpnit) = \J Nn{t) exp{—i(j)n{t)) Can be used as conjugate 
variables with respect to the Hamiltonian iH leading to 
a set of canonical equations idt4'n = dH / dtpn'-i ''■9t4'n = 
-dH/dtpn- Substituting (|4) we get the DNLSE. To sim- 
ulate the leaking process at the two edges, we supplement 
the standard DNLSE with a local dissipation at the two 



edges of the lattice. The resulting leaking DNLSE is [10] 




= (Xl'0n| + /i„)'(/'n - |[V'n-l(l - (^n,!) + V'«+l(l - <5„,Af)] 



- hi'niSn,! + Su^m]; n=l,---,M (5) 

where the normalized time is defined as r = Tt, x = 
2U/T is the rescaled nonlinearity, ^„ = fJ^n/T is the 
rescaled chemical potential and 7 is the atom emission 
probability describing atomic losses. It is useful to de- 
fine the initial effective interaction per well i.e. A = xp 
where p — N{t — 0) /M is the initial average density of 
atoms in the OL. In our numerical experiments we have 
used initial conditions with randomly distributed phases, 
and an almost constant amplitude with only small ran- 
dom fluctuations across the OL. We normalized the wave 
functions such that N{t = 0) = 1. The initial states 
are first thermalized during a conservative (i.e. 7 = 0) 
transient period of, typically, r = 500. The dissipation 
at the lattice boundaries is switched on only after this 
transient is completed, leading to a progressive loss of 
atoms. We study the decay and the statistical properties 
of N{t) as a function of the parameter A. The results 
reported in this Letter correspond to /i„ = and dissi- 
pation rate 7 = 0.2 which is within the experimentally 
accessible range [lO]. Nevertheless we have checked that 
we get the same qualitative behavior for other values of 
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Figure 1: Main panel: Distribution of avalanches 7-'(J) of 
various system sizes M for interatomic interaction strengths 
A — 0.5 and A = 1. In the former case we observe a conver- 
gance to a power law distribution 'P(J) ~ J~" as the lattice 
size M increases, while in the latter case the asymptotic distri- 
bution has already been reached for M = 512. The best least 
square fit indicates that a — 1.86 ± 0.04 in agreement with 
the bounds 1 < a < 3 (see discussion at the text). The distri- 
bution is generated over different initial thermal excitations. 
Upper right inset: Representative realizations of atomic pop- 
ulation decay showing avalanches. Upper left inset: Power 
law distribution of norms P(x = I^Anl^) ~ X-'' for A = 1. The 
best least square fit indicates that f3 — 1.9 ± 0.05 « a. Lower 
left inset: The localization parameter VTZs as a function of 
the initial effective interatomic interaction A. 
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Figure 2: Left panel: Evolution of atomic population for a 
lattice of size M = 128 and various interatomic interaction 
strengths A : (a) A = 0.1 < Ab where no DBs are formed; (b) 
A = 1 > Af, corresponding to the critical regime where scale 
free avalanches are created; (c) A = 16 ^ 1 corresponding 
to the creation of ~ 0{M) DB's (multibreather regime); (d) 
Snapshot of an avalanche event. On the left subpanel, we 
are plotting r vs N{t) whereas on the right we are reporting 
a representative collision event between the outer most self- 
trapped DB and a lattice excitation (moving breather) that 
leads to destabilization of the DB. During the collision, part 
of the moving breather tunnels through the self-trapped DB 
and travels towards the edge of the lattice. The arrival of the 
transmitted density at the edge registers an avalanche in the 
atomic population N{t) (see left subpanel). 
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Figure 3: (a) A destabilization process of a DB hosted by 
a three site system. In the inset we report a Poincare sec- 
tion of the phase space by referring to the variables {N2,'4>), 
where N2 is the norm at site 2 (where the main part of the 
DB is located), and = 1^3 — 02 is the phase missmatch 
between the center, and the next site. In the main figure, 
we report the outgoing atomic population measured at site 
three A^™"^ (associated with an avalanche event-see Fig. 2) 
versus the incoming atomic population Nf^^^ from the first 
site. We observe that atomic population tunnels through 
the DB only if Nf^^ > 0.25, corresponding to the minimal 
excitation needed to trigger the destabilization of the self- 
trapped DB. (b) The distribution 'P(S') of island-sizes for the 
kicked rotor with kicking strength K = 3.5, corresponding to 
a mixed phase space (inset). The island sizes were evaluated 
by studing the evolution of close-by initial conditions as a 
function of their seperation S for succesively longer times r . 
We see that in the limit of r ^ 00 the distribution converge 
to a power law. 



